function [psi,a,ev]=pca(f,N)

% [psi,a,ev]=pca(f,N)
%
% Principal Component Analysis using by the Snapshot Method. (Sirovich, 1987)
% Use SVD instead if (# of images) > (# of pixels/image).
%
% MODEL:  f = a * diag(sqrt(ev)) * psi'
%
% INPUT:
%   f    = matrix of images. (individual images are row vectors)  
%   N    = the number of eigenpairs to compute. [default ALL]
% 
% OUTPUT:
%   psi = column vectors of spatial principal components.
%   a   = comumn vectors of temporal principal components.
%   ev  = PCA eigenvalues in descending order.

if size(f,1) > size(f,2)
  disp('ERROR:  The number of images exceeds the number of pixels/iamge.')
  disp('        Consider using SVD instead.')
  warning off; return;
end

disp('Forming the correlation matrix ...');
cor = f*f';     % the pixel-correlation matrix.

if nargin<2       % If N is not specified, compute all eigenpairs.
  N=size(f,1); 
end;

disp('Diagonalizing ...')
[a,ev]=eig(cor);           % compute the temporal eigenvectors/values.
[ev,ind]=sort(diag(ev));   % sort in ascending magnitude.
ev=flipud(ev);             % switch from ascending to descending.
a=a(:,flipud(ind));        % order the temporal eigenvectors accordingly.
disp('Calculating PCA eigenvectors ...')
a=a(:,1:N);
psi=f'*a*diag(ev(1:N).^-0.5);  % compute the spatial eigenvectors.

